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Abstract. The Laser Interferometer Space Antenna (LISA) is expected to 
detect gravitational radiation from a large number of compact binary systems. 
We present a method by which these signals can be identified and have their 
parameters estimated. Our approach uses Bayesian inference, specifically the 
application of a Markov chain Monte Carlo method. The simulation study that 
we present here considers a large number of sinusoidal signals in noise, and 
our method estimates the number of periodic signals present in the data, the 
parameters for these signals and the noise level. The method is significantly 
better than classical spectral techniques at performing these tasks and does not 
use stopping criteria for estimating the number of signals present. 


PACS numbers: 04.80.Nn, 02.70.Lq, 06.20.Dq 
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1. Introduction 

LISA [Q is expected to detect a very large number of signals from compact binaries 
in the 10 _2 mHz to 100 mHz band, making signal identification very difficult. Tens 
of thousands of signals could be present in the data with significant signal-to-noise 
ratios. In the 0.1 mHz to 3 mHz band there will be numerous signals from white dwarf 
binaries. Sources above 5 mHz should be resolvable, however below ImHz there will 
be source confusion. In the 1 mHz to 5 mHz band we expect as many as 10 5 potential 
sources Guam resulting in an astoundingly difficult data analysis problem. We direct 
the reader to Barack & Cutler [2j and Nelemans et al. [5J for an in-depth description 
of the population of binary systems in the LISA operating band, and how LISA’s 
performance is influenced by them. 

The goal of this paper is to introduce the LISA data analysis community to a new 
approach for identifying and characterizing these numerous signals. We apply Bayesian 
Markov chain Monte Carlo (MCMC) methods to a simplified problem that will serve 
as an example of the technique. MCMC methods are a numerical means of parameter 
estimation, and are especially useful when there are a large number of parameters |BJ. 
We have already applied MCMC methods to other gravitational radiation parameter 
estimation problems; for example, we have used a Metropolis-Hastings (MH) algorithm 
0 El for estimating astrophysical parameters for gravitational wave signals from 
coalescing compact binary systems 0, and pulsars [Him. We believe that MCMC 
methods could provide an effective means for identifying sources in LISA data. We 
summarize our reversible jump MCMC technique in this paper. A more detailed and 
comprehensive description can be found in na 

Here we present a summary of our study of some simple simulated data, 
comprising a number of sinusoidal signals embedded in noise. Our reversible jump 
MCMC algorithm infers the parameters for each (sufficiently large) sinusoidal signal, 
the magnitude of the noise and the number of signals present. In our approach we solve 
both the detection and parameter estimation problems without the need for evaluating 
formal model selection criteria. The method does not require a stopping criterion for 
determining the number of signals and produces results which compare very favorably 
with classical spectral techniques. A Bayesian analysis naturally encompasses Occam’s 
Razor and a preference for a simpler model m- In addition, our MCMC method is 
better than a classical periodogram at resolving signals that are very close in frequency, 
and we provide an explanation of how to identify these signals. 

The method that we present here is not a source subtraction method G3 Signals 
that are sufficiently strong will be identified with a quoted confidence, and sources 
that are weak will simply contribute to the noise, the level of which we also estimate. 
We show that the noise level estimate from our method depends (as it should) on the 
inherent detector noise level, and also the presence of unidentified signals. A benefit 
of using MCMC methods is that computation time does not show an exponential 
increase with the number of parameters 0. 

The problem of identifying an unknown number of sinusoids is neither new nor 
simple [Tin rr> . Previous studies have looked for a handful of unknown signals, here 
we show results for 100 signals. MCMC methods are robust and dynamic, and we 
believe that ultimately it will be possible to use them with LISA data to estimate 
the parameters of all modelled sources types. In the future we will make the model 
more complex, taking into account the orbit of the LISA spacecraft and binary source 
evolution. 
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One can approach the problem of identifying and enumerating sinusoidal signals in 
noise from a number of different directions, but one thing is clear. Discrete noisy 
data can be fitted exactly if one uses a sufficient number of components - the result is 
simply the discrete Fourier transform of the data. Classically, we proceed by estimating 
the noise floor of the spectrum and identify a threshold spectral power that divides 
the components between signals and noise. In this way we prevent the model from 
overfitting the noise. In an iterative fitting procedure, this is achieved by halting when 
the statistics of the residuals fit the noise model well. One very attractive, and well 
known, feature of Bayesian inference is that these ideas are within the fabric of the 
method. Indeed they are such a basic property of logical inference that there is no need 
to refer to ideas such as ‘overfitting’ at all. More generally, the method discourages 
us from using models that have more degrees of freedom than are necessary for the 
problem in hand. 

One can see how this works in a simple example: Take two data models, M\ 
and M 2 constrained by a single datum, d. Mi has one parameter, a\, to describe 
the datum, whereas M 2 uses the sum of two parameters, s = ai + 02 to describe the 
same datum. Which model is better? Here we have no noise and no random variables, 
so this is not a problem for orthodox statistics. However, if the datum is equally 
consistent with both models, we would clearly prefer the simpler model Mi in favour 

of M.2- 

Within the Bayesian framework we consider the odds ratio of the two models: 

_ p(Mi\d) _ p{Mi) p(d\Mi) 

12 p(M2\d) p(M-2) p(d\M.2) 

We will set p(Mi)/p(M 2 ) to unity, as we have no prior preference for either model, 
and take the priors for a\ and <22 to be each uniform in the range 0 —> R , making the 
prior for s the convolution of two of these. The functional priors for a\ under Mi and 
s under M .2 are therefore 


p(ai\Mi) 


1 _ 
R ’ 


p(s\M 2 ) 


s/R 2 for 0 < s < R 

2/R-s/R 2 for R < s < 2 R. 


( 2 ) 


The probability of the data, given either model, is simply a Dirac delta function 
centred on the value of the datum, so we can calculate the evidences p(d\Mi) and 
p(d\M 2 ) by marginalising over the allowed parameter values: 


p(d\Mi) 

p(d\M 2 ) 


J S(d - ai) dai = j 

r2R 

/ p(s\M 2 )S(d — s) ds = 

Jo 


1/R 

0 


for 0 < d < R 
otherwise, 

d/R 2 for 0 < d < R 

2/R - d/R 2 for R<d<2R 
0 otherwise. 


( 3 ) 

( 4 ) 


as shown in Fig. |T] If the datum lies in the range 0 < d < R, the odds ratio is 
(1/R)/(d/R 2 ) = R/d , so that Mi will always favoured over M 2 in that range. If 
d = R the odds ratio is unity, and if R < d < 2R the only possible model is M 2 - 

This demonstrates why the simpler model is favoured even when the datum is 
equally consistent with both: M 2 has more flexibility than is necessary to explain the 
datum and so penalizes itself by spreading its evidence more thinly. 
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Figure 1. The evidences for models Adi (solid line, one parameter fit) and M 2 
(dashed line, two parameter fit) as a function of the datum d. Note that the 
evidence ratio always favours the simpler model Mi when the datum is equally 
consistent with either. 


Although in this example we consider a probability ratio to determine our 
favoured model, when more than two models are available we can consider the model 
choice to be a parameter itself, and determine its marginal posterior probability in 
the usual way. 

A slightly more pertinent, though still highly restricted, example would be a 
data set {d k } that consists of observations of the sum of m sinusoids of the form 
Ai sin(27r/)t) at times t = t k and Gaussian noise with variance tJ 2 . We also know 
that 0 < to < 5, and Ai and fi can only take the discrete values A,; £ {1,..., 5} and 
fi £ {0.01,0.02, ...,0.05}. 

For this discrete problem, assuming uniform priors on m, Ai and fi for each 
sinusoid, we can identify the probability of any particular to, given the data, 
irrespective of the other signal parameters: 

p(m\{d k }) oc ^ H ex P(-X 2 / 2 ) (5) 

j4i,/i AmJm 


where 


X 




k 


I 2 


d k - ^ A* sin(27r/,4) 


( 6 ) 


Here it is the factor of 25 m , originating from the to normalised priors for the model 
parameters, that offsets the increasingly good y 2 fit that might come from large values 
of to and provides our Occam factor. 

This discrete problem can be solved by directly marginalising over the nuisance 
amplitude and frequency parameters. However problems with more and/or with 
continuous parameters are not approachable using such direct methods, and must 
be tackled in another way. 


3. Parameter estimation 

We consider the continuous case as a signal consisting of to superimposed sinusoids 
where m is unknown. Therefore we confine our attention to a set of models 
{M m : to £ {1, • • •, M}} where M is the maximum allowable number of sinusoids. 
Let d = be a vector of N samples recorded at times t = [ii, • • • 

Model M m assumes that the observed data are composed of a signal plus noise: 


1 
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dj = a rn ) + 6j , for j = 1... N , where the noise terms ej are assumed to be i.i.d. 
N( 0, a^f) random variables. The signal of model A4 m is assumed to be of the form 


m 

f m(tj,a m ) = J2 [ A i m) cos(2t rf z (m) tj) + B f m) sin(27r/| m) t j ) 


Model M m is therefore characterized by a vector 

n — f( m ) ... /l( m ) R( m ) f( m ) n- 2 1 

“m “ 1^1 J- L> 1 )J1 ) ) iJm 


(7) 

( 8 ) 


of 3m + 1 unknown parameters. The objective is to find the model A4 m that best fits 
the data. To this end, we use a Bayesian approach as in El The joint probability of 
these data d given the parameter vector a m and model A4 m is 


p(d\m, a m ) oc — w exp 

a m 


1 

2 ^m 


N 2 | 

'y ' [dj ~ f m{tj, dm)} / 

i = 1 J 


(9) 


We choose (now continuous) uniform priors for the amplitudes and 

frequency f\ m) with ranges A[ m) e [-A max ,^ max ], £ [-S max ,B max ] and 

^ £ [0, 0.5], respectively. Furthermore we use a uniform prior for m over {1,..., M} 
and vague inverse Gamma priors for By applying Bayes’ theorem, we obtain the 
posterior pdf 


p(m, a m \d) 


p{m, a m )p{d\m, a m ) 
p(d) 


( 10 ) 


where p(d) = IP( m i a m)p{d\m, a m ) d a m . We use a sampling-based technique 

for posterior inference via MCMC 'Bjj. MCMC techniques only require the 
unnormalised posterior p(m, a m \d) oc p(m , a m )p(d\m, a m ) to simulate from Eq. (11011 
in order to estimate the quantities of interest. However, as the dynamic variable of the 
simulation does not have fixed dimension, the classical MH techniques 0IHI cannot be 
adopted when proposing trans-dimensional moves between models where the model 
indicator m determines the dimension (3m + 1) of the parameter vector a m . We 
therefore use the Reversible Jump Markov Chain Monte Carlo (RJMCMC) algorithm 
CHI IE] for model determination, as in m • For transitions within the same model, 
we use the delayed rejection method EHIE] which yields a better adaptation of the 
proposals in different parts of the state space. 


3.1. The RJMCMC for model determination 

To sample from the joint posterior p(rn,a m \d) via MCMC, we need to construct a 
Markov chain simulation with state space U^f =1 (m x JR 3m+1 ). When a new model is 
proposed we attempt a step between state spaces of different dimensionality. Suppose 
that at the nth iteration of the Markov chain we are in state ( k,a,k ). If model Mk 1 
with parameter vector a' k , is proposed, a reversible move has to be considered in 
order to preserve the detailed balance equations of the Markov chain. Therefore 
the dimensions of the models have to be matched by involving a random vector r 
sampled from a proposal distribution with pdf q , say, for proposing the new parameters 
a' k , = t(ofc, r) where t is a suitable deterministic function of the current state and the 
random numbers. Here we focus on transitions that either decrease or increase models 
by one signal, i.e. k' £ {k — l,k + 1}. We use equal probabilities Pk^k’ = Pk’^k to 
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either move up or down in dimensionality. Without loss of generality, we consider 
k < k!. 

If the transformation t k \^k' from (a k ,r) to a' k , and its inverse thu = t k’t-^k 
are both differentiable, then reversibility is guaranteed by defining the acceptance 
probability for increasing a model by one signal according to m by 


a k ^k'{a' k ,\cL k ) = min 


l p(a' k ,, k')p(d\a' k ,, k')p k ^ k , 

’ p{a k , k)p(d\a k , k)q(r)p k ^ k 


\Jk~k’\ (11) 


where \Ji 


k\-^k' | — 


is the Jacobian determinant of this transformation and 


dt(a k ,r) 
d(a k ,r) 

q(r) is the proposal distribution. 

In this context, two types of transformations, ‘split-and-merge’ and ‘birth-and- 
death’, are obvious choices. In a ‘split-and-merge’ transition, the proposed parameter 
vector a' k , comprises all (k — 1) subvectors of a k except a randomly chosen subvector 

(k~) (b) (k) 

a = (A\ , B\ , fl 'y which is replaced by two 3-dimensional subvectors, 


°(u) = 



p(k ) 
’Wi 



and 


“fe) - 



■ II 


(k’) 



with roughly half the amplitudes but about the same frequency as ap). 

A three-dimensional Gaussian random vector (with mean zero), r = (rA,rB,rf), 
changes the current state a^ to the two resulting states a'^, through a linear 
transformation 


t fc^fc' (a (i) , r) = + r A , } + r B ,f[ k) + ry, ^A [k) - r A ,\B {k) - r Bl f [k) - r f ^j 

_ ( 4(^0 0(^0 a(^') ll(k') r(k')\ 

’ ^ii ’ Jii ’ J12 J ‘ 


( 12 ) 


By analogy, the inverse transformation tjWw := t k '<-, k accounts for the merger of two 
signals. Note that the determinant of the Jacobian of the transformation t k ^, k > is 
| = 2, and that of its inverse is 1/2. 

We use a multivariate normal distribution, N[0, diag(cr^, a^, cr 2 )], for the proposal 
distribution q(r). Care has to be taken in choosing suitable values for the proposal 
variances to achieve reasonable acceptance probabilities in Eq. CU). 

The second ‘birth-and-death’ transformation consists of the creation of a new 
signal with parameter triple a/- independent of other existing signals in the current 
model A4 k . The one-to-one transformation in this case is very simply given by 
t kk ^ k '(r) = r = ct( fc+1 ) with Jacobian equal to 1. 

Here, q is the three-dimensional pdf from which we draw proposals for the 
additional signal. We use independent uniform distributions with frequency range 
0 < / < 0.5 and amplitude range ( Af + B 2 ) 1 / 2 < l a where l a is the radius for the two 
amplitudes of the signal in polar coordinates. Again, the radii of the uniform proposal 
densities have to be tuned to achieve an optimal acceptance rate. 


3.2. The delayed rejection method for parameter estimation 

For transitions within a model A4 m , classical MCMC methods can be applied. 
Here, however, we use an adaptive MCMC technique, the delayed rejection method 
EDI EH EH that we have successfully applied to estimate parameters of pulsars EH 
The idea behind the delayed rejection method is that persistent rejection indicates 
that locally, the proposal distribution is badly calibrated to the target. Therefore, 






LISA source confusion: identification and characterization of signals 


7 


the MH algorithm is modified so that on rejection, a second attempt to move is made 
with a proposal distribution that depends on the previously rejected state. In this 
context, when a proposed MH move is rejected from a bold normal distribution with 
large variance a second candidate can be proposed with a timid proposal distribution 
for sampling the parameters for the individual sinusoids. Hence, the main objective 
of the first stage is a coarse scan of the parameter space and therefore we choose 
the variances of the parameters about one order of magnitude smaller than the prior 
ranges of the corresponding parameters. Once a mode is found, we aim to draw 
representative samples in the second stage. 

The precision of the frequency in a single-frequency model depends on the 
amplitude, the variance a 1 of the noise, and the number of samples N of the data 
set [321122. The precision of the frequency has been derived in m by a Gaussian 
approximation to the posterior pdf of the frequency and calculation of its standard 
deviation, given by a'j = (27r) _1 [48cr 2 /7V 3 (H 2 -|-I3 2 )] 1 / 2 . We therefore choose proposals 
with this standard deviation. 

3.3. Starting values 

The starting values of a Markov chain are crucial for the length of the burn-in 
period, i.e. the time needed for the chain to achieve convergence to the real posterior 
distribution. We perform a Fast Fourier Transformation (FFT) prior to the simulation 
and use corresponding estimates as starting values. Arthur Schuster introduced the 
periodogram m 

C(f) = ^ [R(f) 2 + 1(f) 2 ] (13) 

where R(f) = Y^jLi dj cos (27r/tj) and 1(f) = dj sin (2irftj) are the real 

and imaginary parts from the sums of the discrete Fourier transformation. As a 
starting value for m, we use the number mo of local maxima in the periodogram that 
exceed a certain noise level (lower than the expected one). We use the frequencies 
corresponding to the local maxima in the periodogram, /oy, as starting values for 
/q 7 o) and A 0 y = 2R(fo l i)/N, B 0ti = 2I(f 0 ^)/N as starting values for A- m °^ and 
B[ m °\ respectively. 

3.fi Identifying the sinusoids 

Although the RJMCMC offers great ease in model selection, we still encounter the 
label-switching problem, a general problem due to invariance of the likelihood under 
relabelling that has been extensively discussed in the context of mixture models Eg 

The sinusoids that are contained in the model m with the highest posterior 
probability of m are permutations of m coexistent sinusoids out of a number of 
sinusoids that we do not know but can estimate by the upper limit m max of the 
marginal posterior of m. Therefore, the parameter vector sampled in each iteration of 
the Markov chain (corresponding to model m) is a permutation of m parameter triples 
determining to out of m max sinusoids. The problem is to determine which parameter 
triple belongs to which sinusoid. 

The parameter that contributes significantly to identifying a sinusoid is its 
frequency. We thus calculate the marginal posterior of the frequency and obtain the 
mmax strongest peaks together with their frequency ranges by finding the threshold 
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that separates those peaks. It is still possible that individual peaks contain more than 
one sinusoid or even none. This can be assessed by a histogram simular to that in Fig. [21 
but restricted to the frequency range under consideration. To separate more than one 
present sinusoids, we then consider the two amplitudes and apply an agglomerative 
hierarchical cluster analysis that involves all three parameters. We use a modified 
Ward technique m that minimizes the within group variance using a normalized 
Euclidean distance between the parameters by adjusting the frequency range to the 
much larger range of the amplitudes. 

4. Results 

We created an artificial data set of 1000 samples from m = 100 sinusoids. The 
sinusoids were randomly chosen with maximum amplitudes 1 and the noise standard 
deviation was a = 1. We chose a uniform prior for m on {0, 1,2,... ,M = 60000}, 
and set A max = R max = 5. The Markov chain ran for 10 8 iterations and was thinned 
by storing every 1 000th iteration. The first 5 x 10 6 iterations where considered as 
burn-in and discarded. The MCMC simulation was implemented in C on a 2.8 GHz 
Intel P4 PC and took about 43 hours to run. Fig. [2] gives the histogram of the 
posterior model probabilities obtained by the reversible jump algorithm. As each 
model M m is characterized by a different noise level er m , we have also plotted the 
posterior distributions of the noise standard deviations for increasing model order. 
Note that rr m decreases with higher model order m since a model comprising more 


Posterior model probabilities Median and 95%-confidence intervals of a for different models 



Figure 2. Posterior distribution of model number m with 5%, 50%, and 95% 
posterior confidence intervals of the corresponding noise levels. 

sinusoids accounts for more noise. Here, we choose model m = 95 corresponding to 
the posterior mode of m as the best fitting model. 

We used all MCMC samples corresponding to model Mg 5 . For 

ease of notation, we denote the parameter vector of model .M 95 by 
(A 1; Bi, fi,..., A 95 , B g5 , / 95 ). The complete line power spectrum density can 
be estimated by the product of the conditional expectation of the energy 
E{A\ + Bf\d, m, fi) of each sinusoid * given its frequency fi, and the posterior pdf 
of fi given the data, p(fi\d). One of the advantages of the Bayesian spectrum anal¬ 
ysis is the possibility to calculate confidence areas for the spectrum. Therefore we 
group our MCMC samples and calculate posterior confidence intervals for each fre¬ 
quency bin. A sufficient width for the bins can be assessed by the frequency accuracy 
af = (27r • snr) -1 (48/IV 3 ) 1 ' /2 given by {T7j , where ‘snr’ is the signal-to-noise ratio. In 
our example, the choice of 30 000 bins is sufficient to resolve sinusoids with an snr of 
about 2. Fig. [21 shows the real signals, the Bayes spectrum and the classical Schuster 
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Real signals 




Periodogram 



Figure 3. Comparison of real signals, classical Schuster periodogram, and 
Bayesian spectrum. Note that the total energy of the three spectra is similar 
but the lines have different accuracies. Therefore the energy of a single line is 
more concentrated due to the better accuracy of the Bayesian spectrum estimate 
resulting in higher peaks. The theoretical spectrum consists of delta functions for 
each sinusoid which would yield infinite peaks. Therefore, an energy contribution 
of ^ (A? + B ?) for each sinusoid with frequency fi is shown. 

The plot for the real signals displays an individual energy contribution for each 
sinusoid i of (A? + Bf)N/2. Normally a theoretical spectrum would consist of 
delta functions with infinitely large energy peaks since the energy contribution is 
concentrated on an interval of infinitely small width. Therefore we just plotted the 
energy contribution on the energy scale that yields a similar scaling as obtained by 
the periodogram. 

In order to be able to display 95% confidence areas of the spectrum we present 
three magnified areas in Fig.0J In plot (a) we see a sinusoid with rather small energy. 
The accuracy of the frequency estimation is worse compared with the sinusoid of 
graph (b) which has a significantly larger energy. The third graph shows two very 
close sinusoids. The frequency estimation is very inaccurate due to the interference of 
the two signals. This is consistent with theoretical results by E The interference of 
the two close signals is due to a phase shift of 175°. The interference and hence the 
frequency estimation depends upon whether the sinusoids are orthogonal E or not. 
Nevertheless, we are able to identify the existence of two signals while the periodogram 
only reveals the existence of one. 
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(a) Median and 95% confidence 
area of energy for sinusoid: 
E=78.2 [45.96,119.9] 

Energies from corresponding 
real parameters: 

E=49.9 at f=0.0500462 



(b) Median and 95% confidence 
area of energy for sinusoid: 
E=476.1 [389.9,570.7] 

Energies from corresponding 
real parameters: 

E=459.5 at f=0.162425 



(c) Median and 95% confidence 
area of energy for sinusoids: 
E=3497.2 [819.7,10800.8] 
E=2731.6 [261.8,8809.1] 
Energies from corresponding 
real parameters: 

E=462.7 at f=0.360133 
E=66.9 at f=0.36065 


Figure 4. Three magnified areas for a comparison of the classical periodogram 
(dotted line) and the 95% confidence area (area within bright colored lines) of the 
Bayesian spectrum. The darker colors within the 95% confidence area indicate the 
median of the energy. The vertical lines show the delta functions of the theoretical 
spectrum. Its energies are given above. 


The estimates of the amplitudes, however, always show huge values and confidence 
intervals for sinusoids close in frequency. The huge energies are merely restricted by 
the choice of priors for the amplitudes. The reason for this is due to the possible 
combinations to express a sum of sinusoids when the observation time is insufficiently 
long with respect to the distance in frequency. In this case we can not make accurate 
statements about the amplitudes and hence energies of both sinusoids. 

If we take a look at the single peak of the periodogram the energy that is 
considered is subject to the data from a discrete and finite observation time, given 
by This, however does not reflect the energy contribution of the 

real signals. The Bayesian estimates of the amplitudes are honest by yielding large 
confidence areas for the energies of sinusoids close in frequencies but in return small 
confidence intervals for isolated sinusoids. 

5. Discussion 

We have presented a Bayesian approach to identifying a large number of unknown 
periodic signals embedded in noisy data. A reversible jump MCMC technique can be 
used to estimate the number of signals present in the data, their parameters, and the 
noise level. This approach allows for simultaneous detection and parameter estimation, 
and does not require a stopping criterion for determining the number of signals. The 
MCMC method compares favorably with classical spectral techniques. 

Our motivation for this research is to address the difficulty that LISA will 
ultimately encounter in having too many signals present. LISA may see 100 000 signals 
from binary systems in the 1 mHz to 5 mHz band. We see our work as a new method 
that could help LISA to identify and characterize these signals. The work here is a 
simplified problem, one that neglects the time evolution of the signal and modulation 
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due to LISA’s orbit. The next step is to deal with these more complicated signals, and 
to develop a realistic strategy for applying our MCMC methods to more realistic LISA 
data. We believe that MCMC methods, like those presented here, provide a practical 
and highly effective method of identifying and characterizing the large number of 
signals that will exist in the LISA data. 
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